Plotting 39 Storms’ Precipitation for First 24 Hours After Landfall

all.storm.folders <- list.dirs("NAMandST4", recursive = F) #G:\\NAMandST4
NAMlist <- list()
ST4list <- list()
mask <- raster("G:\\lsmask.nc")
mask[mask==-1] <- 0
extent(mask)[1] <- extent(mask)[1]-360
extent(mask)[2] <- extent(mask)[2]-360
mask.regrid <- resample(mask, projectRaster(raster("G:\\nam_218_20050829_1200_f012.grib"),crs = "+proj=longlat +datum=WGS84"), method='ngb')


for(i in 1:length(all.storm.folders)){
  # [-c(15:36)],[c(1:16,35:39)] right after all.storm.folders to use only the ones without jpeg2000 issue
  # i <- 1
  print(i)
  currentST4_NAM_folders <- list.dirs(all.storm.folders[i])[-1] #uncomment the [-1] when st4 and nam folders are together
  ST4_folder <- currentST4_NAM_folders[1]
  NAM_folder <- currentST4_NAM_folders[2] #uncomment the [2] as well
  
  print(length(list.files(NAM_folder)))
  #print(list.files(NAM_folder))
  NAMlist[[i]] <- list()
  ST4list[[i]] <- list()
  if(length(list.files(NAM_folder,full.names = T))==18){first24 <- c(7,8)} else 
  if(length(list.files(NAM_folder,full.names = T))==15){first24 <- c(5,6)} else {print("something went wrong"); next}
  for (j in 1:length(list.files(NAM_folder,full.names = T)[first24])) {
    # j <-1
    NAMfile <- list.files(NAM_folder,full.names = T)[j]
    # paste(NAM_folder,"/",NAMfile, sep = "") #don't need when full.names=T in list.dirs
    date_start <- str_locate_all(pattern ='_218_', NAMfile)[[1]][1,2] + 1
    storm_date <- as.numeric(substr(NAMfile,date_start,date_start+7))
    if(storm_date < 20060701){NAM_band <- 10} else if(storm_date > 20060701 && storm_date < 20170101){NAM_band <- 11} else {NAM_band <- 373}
    print(NAM_band)
    NAMlist[[i]][[j]] <- (raster(NAMfile, band=NAM_band) %>% projectRaster(crs = "+proj=longlat +datum=WGS84"))
    print(min(values(NAMlist[[i]][[j]]), na.rm=T))
    values(NAMlist[[i]][[j]])[which(values(NAMlist[[i]][[j]]<0))] <- 0 #change negative projected precipitation to 0
  }
  for (k in 1:length(list.files(ST4_folder,full.names = T)[1:4])) {
    # k <-1
    ST4file <- list.files(ST4_folder,full.names = T)[k]
    ST4list[[i]][[k]] <- (raster(ST4file) %>% projectRaster(crs = "+proj=longlat +datum=WGS84"))
    ST4list[[i]][[k]] <- resample(ST4list[[i]][[k]], mask.regrid, method = "bilinear")
    
    print(min(values(ST4list[[i]][[k]]), na.rm=T))
    values(ST4list[[i]][[k]])[which(values(ST4list[[i]][[k]]<0))] <- 0 #change negative projected precipitation to 0
  }
  first24NAM <- (NAMlist[[i]][[1]]+NAMlist[[i]][[2]])*mask.regrid
  first24ST4 <- (ST4list[[i]][[1]]+ST4list[[i]][[2]]+ST4list[[i]][[3]]+ST4list[[i]][[4]])*mask.regrid
  error.max <- max(abs(min(values(first24NAM-first24ST4),
                           na.rm = T)),max(values(first24NAM-first24ST4),na.rm = T))

  plot(first24NAM, xlim=c(-110,-55), ylim=c(23,50),
       main=paste(all.storm.folders[i],"NAM"));plot(st_geometry(spData::us_states), add=TRUE)
  plot(first24ST4, xlim=c(-110,-55), ylim=c(23,50),
       main=all.storm.folders[i]);plot(st_geometry(spData::us_states), add=TRUE)
  plot(first24NAM-first24ST4, breaks=seq(-1*ceiling(error.max/10)*10,ceiling(error.max/10)*10,ceiling(error.max/10)),
       xlim=c(-110,-55), ylim=c(23,50),col=colorspace::diverge_hsv(20),
       main="Error, Forecast NAM - Observed ST4")
  plot(st_geometry(spData::us_states), add=TRUE)
}
## [1] 1
## [1] 15
## [1] 10
## [1] -0.9064407
## [1] 10
## [1] -0.270989
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 2
## [1] 15
## [1] 10
## [1] -0.4068737
## [1] 10
## [1] -0.3602461
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 3
## [1] 15
## [1] 10
## [1] -0.2420354
## [1] 10
## [1] -0.5602488
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 4
## [1] 15
## [1] 10
## [1] -0.3504511
## [1] 10
## [1] -0.3234572
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 5
## [1] 15
## [1] 10
## [1] -0.2141921
## [1] 10
## [1] -0.5485855
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 6
## [1] 15
## [1] 10
## [1] -0.4267543
## [1] 10
## [1] -0.3410163
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 7
## [1] 15
## [1] 10
## [1] -0.8998696
## [1] 10
## [1] -1.492411
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 8
## [1] 15
## [1] 10
## [1] -0.1717761
## [1] 10
## [1] -0.6953588
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 9
## [1] 15
## [1] 10
## [1] -0.1503301
## [1] 10
## [1] -0.09110766
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 10
## [1] 15
## [1] 10
## [1] -0.3714533
## [1] 10
## [1] -0.1491003
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 11
## [1] 15
## [1] 10
## [1] -0.1790766
## [1] 10
## [1] -0.451881
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 12
## [1] 15
## [1] 10
## [1] -0.04831516
## [1] 10
## [1] -0.3113824
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 13
## [1] 15
## [1] 11
## [1] -0.6925985
## [1] 11
## [1] -0.6240735
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 14
## [1] 15
## [1] 11
## [1] -0.7908002
## [1] 11
## [1] -0.680807
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 15
## [1] 18
## [1] 11
## [1] -0.3895373
## [1] 11
## [1] -0.1227917
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 16
## [1] 18
## [1] 11
## [1] -0.442489
## [1] 11
## [1] -0.3708847
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 17
## [1] 18
## [1] 11
## [1] -0.4058006
## [1] 11
## [1] -0.5956256
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 18
## [1] 18
## [1] 11
## [1] -0.4088998
## [1] 11
## [1] -0.4297405
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 19
## [1] 18
## [1] 11
## [1] -0.4451192
## [1] 11
## [1] -0.1183022
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 20
## [1] 18
## [1] 11
## [1] -0.5069734
## [1] 11
## [1] -0.7084994
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 21
## [1] 18
## [1] 11
## [1] -0.3589332
## [1] 11
## [1] -0.1560499
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 22
## [1] 18
## [1] 11
## [1] -0.4813572
## [1] 11
## [1] -0.2920385
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 23
## [1] 18
## [1] 11
## [1] -0.2158466
## [1] 11
## [1] -0.3488085
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 24
## [1] 18
## [1] 11
## [1] -0.282494
## [1] 11
## [1] -1.330304
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 25
## [1] 18
## [1] 11
## [1] -0.9493138
## [1] 11
## [1] -1.91404
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 26
## [1] 18
## [1] 11
## [1] -0.1666827
## [1] 11
## [1] -0.09597836
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 27
## [1] 18
## [1] 11
## [1] -0.3948196
## [1] 11
## [1] -0.2775595
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 28
## [1] 18
## [1] 11
## [1] -0.9548955
## [1] 11
## [1] -0.06103402
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 29
## [1] 18
## [1] 11
## [1] -0.1150618
## [1] 11
## [1] -0.05929433
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 30
## [1] 18
## [1] 11
## [1] -0.3006031
## [1] 11
## [1] -0.1595173
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 31
## [1] 18
## [1] 11
## [1] -0.3618694
## [1] 11
## [1] -0.286569
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 32
## [1] 18
## [1] 11
## [1] -0.06227638
## [1] 11
## [1] -0.1600004
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 33
## [1] 18
## [1] 11
## [1] -0.1771795
## [1] 11
## [1] -0.1172347
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 34
## [1] 18
## [1] 11
## [1] -0.1443691
## [1] 11
## [1] -0.6040375
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 35
## [1] 18
## [1] 11
## [1] -1.115945
## [1] 11
## [1] -0.06179695
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 36
## [1] 18
## [1] 11
## [1] -0.1229631
## [1] 11
## [1] -0.0505522
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 37
## [1] 18
## [1] 373
## [1] -0.3943033
## [1] 373
## [1] -0.5884688
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 38
## [1] 18
## [1] 373
## [1] -0.1104949
## [1] 373
## [1] -2.127217
## [1] 0
## [1] 0
## [1] 0
## [1] 0

## [1] 39
## [1] 18
## [1] 373
## [1] -0.2599918
## [1] 373
## [1] -0.1760313
## [1] 0
## [1] 0
## [1] 0
## [1] 0